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INTRODUCTION 


There  Is  an  ever  Increasing  need  to  solve  problems  of  greater  complexity 
and  a  corresponding  need  for  reliable  and  robust  software  tools  to  accurately 
and  efficiently  describe  the  phenomena.  Adaptive  techniques  are  good 
candidates  for  providing  the  computational  methods  and  codes  necessary  to 
solve  some  of  these  difficult  problems.  Two  popular  adaptive  techniques  are: 
(1)  moving  mesh  methods,  where  a  grid  of  a  fixed  number  of  finite  difference 
cells  or  finite  elements  is  moved  in  order  to  follow  and  resolve  local 
nonuniformities  in  the  solution,  and  (2)  local  refinement  methods,  where 
uniform  fine  grids  are  added  to  coarser  grids  in  regions  where  the  solution  is 
not  adequately  resolved.  A  representative  sample  of  both  types  of  mechods  is 
contained  in  Babuska,  Chandra,  and  Flaherty  (ref  1).  Recently,  Adjerid  and 
Flaherty  (ref  2)  developed  a  finite  element  method  that  combines  mesh  moving 
and  refinement. 

Herein,  we  discuss  a  local  refinement  finite  element  procedure  for 
finding  numerical  solutions  of  M-dimensional  vector  systems  of  partial 
differential  equations  having  the  form 

Lu  ut  +  f(x,t,u,ux)  -  [D(x,t ,u,ux]x  -0  ,  a<x<b,  t>0  (1) 

subject  to  the  initial  conditions 

u(x,0)  ■  u°(x)  ,  a  <  x  <  b  (2) 

and  appropriate  boundary  conditions  so  that  the  problem  has  a  well-posed 
solution. 


^-Babuska,  I.,  Chandra,  J.,  and  Flaherty,  J.  E.  (eds.).  Adaptive  Computational 
Methods  for  Partial  Differential  Equations,  SIAM,  Philadelphia,  1983. 
^Adjerid,  S.  and  Flaherty,  J.  E.,  "A  Moving  Finite  Element  Method  for  Time 
Dependent  Partial  Differential  Equations  with  Error  Estimation  and 
Refinement,"  to  appear  in  SIAM  J.  Numer.  Anal.,  1985. 


We  discretize  Eqs.  (1)  and  (2)  for  a  time  step  using  a  finite  eleraent- 
Galerkin  procedure  with  piecewise  bilinear  approximations  on  a  rectangular 
space-time  net.  At  the  end  of  each  time  step,  we  estimate  the  local 
discretization  error,  add  finer  subgrids  of  space-time  elements  in  regions  of 
high  error,  and  recursively  solve  the  problem  again  in  these  regions.  The 
process  terminates  when  the  error  estimate  on  each  grid  is  less  than  a 
prescribed  tolerance.  The  original  coarse  space-time  grid  is  then  carried 
forward  for  the  next  time  step  and  the  strategy  is  repeated.  Our  algorithm  is 
discussed  further  in  Flaherty  and  Moore  (ref  3)  and  some  of  this  discussion  is 
repeated  in  Che  following  section. 

Berger  (ref  4)  used  a  similar  local  refinement  procedure  to  solve  one- 
and  two-dimensional  hyperbolic  systems.  She  used  explicit  finite  difference 
schemes  to  discretize  the  partial  differential  equations,  while  we  use 
implicit  finite  element  techniques  since  we  are  primarily  interested  in 
parabolic  problems. 

In  addition  to  the  discretization  technique,  the  major  numerical 
questions  that  must  be  answered  as  part  of  the  development  of  a  local 
refinement  code  are  (1)  the  estimation  of  the  discretization  error,  and  (2) 
the  appropriate  initial  and  boundary  conditions  to  apply  at  coarse-fine  mesh 
interfaces.  Of  course  computer  science  questions,  such  as  which  language  to 

^Flaherty,  J.  E.  and  Moore,  P.  K. ,  "An  Adaptive  Local  Refinement  Finite 
Element  Method  for  Parabolic  Partial  Differential  Equations,”  Proceedings  of 
the  International  Conference  on  Accuracy  and  Estimates  and  Adaptive 
Refinements  in  Finite  Element  Computations,  Technical  University  of  Lisbon, 
Lisbon,  Vol.  2,  1984,  pp.  139-152. 

^Berger,  M.  J.,  "Adaptive  Mesh  Refinement  for  Hyperbolic  Partial  Differential 
Equations,”  Report  No.  STAN-CS-82-924,  Department  of  Computer  Science, 
Stanford  University,  1982. 
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use  to  describe  and  Implement  the  various  algorithms  and  what  data  structures 
to  use  to  represent  and  store  the  grids  and  solutions,  must  also  be  answered. 
Our  work  in  all  of  these  areas  is  still  far  from  complete,  and  herein  we  only 
discuss  our  progress  and  thoughts  on  error  estimation  techniques,  data 
structures,  and  interface  conditions  (see  the  following  section).  We  also 
present  the  results  of  three  examples  that  illustrate  our  method  and  some 
preliminary  conclusions  and  future  plans. 


FINITE  ELEMENT  ALGORITHM 

We  discretize  Eq.  (1)  on  a  strip  ct  <  x  <  0,  p<t<q  using  a  finite 
element-Galerkin  method  with  a  uniform  grid  of  N  rectangular  elements  of  size 
(0-a)/N  by  (q-p).  We  refer  to  this  grid  as  R(  a, 0,p,q,N,f ,s) ,  where  f  and  s 
are  pointers  to  the  father  and  son  grids  discussed  later.  Each  grid  uses 
records  to  store  the  appropriate  information. 

We  generate  the  discrete  system  on  R( a, 0,p,q,N,f ,s)  in  the  usual  manner; 
thus,  we  approximate  u  by  U(x,t)  and  select  test  functions  V(x,t),  where  U  and 
V  are  elements  of  a  space  of  C°  bilinear  polynomials  with  respect  to  the  grid 
R.  We  then  take  the  inner  produce  of  Eq.  (1)  and  V,  replace  u  by  U,  and 
integrate  any  diffusive  terms  by  parts  to  obtain 

J  [VTUt  +  vTf(x,t,U,Ux)  +  V^D(x,t ,U)UX] dxdt 
R  x 


vTo(x,t,U)Ux|  dt 


0 


(3) 


Equation  (3)  must  vanish  for  all  bilinear  functions  V  on  the  grid  R.  The 
integrals  are  approximated  using  a  four-point  Gauss  quadrature  rule  and  the 
resulting  nonlinear  system  is  solved  by  Newton  iteration  (see,  for  example 


Reference  5  for  additional  details) .  Appropriate  initial  and  boundary 
conditions  for  Eq.  (3)  are  discussed  later  in  this  section. 

We  describe  our  local  refinement  procedure  for  solving  Eqs.  (1)  and  (2) 
for  one  time  step  (t^t1)  on  a  coarse  grid  with  N°  elements,  i.e.,  on 
R(a(b,t°(t1tN°,Ots)  (where  the  pointer  f  -  0  signifies  that  this  grid  has  no 
father).  To  solve  this  problem,  we  simply  call  the  procedure  "locref"  with 
the  arguments  R(a,b,t°,t *,N°,0,s) ,  tol,  tsub  for  each  coarse  grid  time 
interval.  A  pseudo-PASCAL  description  of  the  procedure  "locref”  is  shown  in 
Figure  1. 


PROCEDURE  locref  (R(a, 3,p,q,N,f ,s) ,  tol,  tsub) 

BEGIN 

Solve  the  finite  element  equations  (Eq.  (3))  on  R(  a,  3,p,q,N,f ,s); 
Estimate  the  error  on  R(  a,  S,p,q,N,f ,s) ; 

IF  error  >  tol,  THEN 
BEGIN 

calculate  where  error  >  tol  and  return  the  son  grids; 

FOR  j  1  TO  tsub  DO 

FOR  i  j-  1  TO  number  of  sons  DO 
BEGIN 

P[JJ  P  +  ( j-i)*(q-p)/tsub; 

HU]  Plj]  +  (q~p)/ tsub; 

locref  (R( <*[i] ,  S[i]  ,p[  j]  ,q[  j]  ,N[i] , 

R(a»  B,p,q,N,f ,8) ,s[i] , tol, tsub) 

END 

END 

END; 


Figure  1.  Algorithm  for  local  refinement  solution  of  Eqs.  (1)  and  (2)  on 

R(“»3»P,q,N,f ,s)  with  an  error  tolerance  of  tol  and  dividing  the 
local  time  step  by  tsub  each  time  the  error  test  is  not 
satisfied. 


The  recursive  algorithm  locref  sets  up  a  tree  structure  of  grids  with 
R(a,b,t°,t 1 ,N° ,0,s)  being  the  root  node  and  with  the  solution  being  generated 


Davis,  S.  F.  and  Flaherty,  J.  E.,  "An  Adaptive  Finite  Element  for  Initial- 
Boundary  Value  Problems  for  Partial  Differential  Equations,"  SIAM  J.  Sci. 
Stat.  Comout..  Vol.  3,  1982,  pp.  6-27.  - 


4 


by  a  preorder  traversal  of  the  tree  at  each  local  time  step.  For  example,  if 
the  root  grid  is  refined  to  given  two  subgrids  and  the  time  step  is  halved, 
the  problem  is  solved  on  the  first  subgrid  on  its  first  time  step,  followed  by 
the  second  subgrid  on  the  same  time  step.  This  procedure  is  repeated  for  the 
second  time  step.  The  error  is  estimated  by  Richardson  extrapolation,  i.e., 
the  space  and  time  steps  are  halved  and  the  problem  is  solved  again  on  this 
new  grid.  The  two  solutions  that  are  obtained  at  each  original  grid  point  are 
used  to  generate  an  error  estimate.  If  this  pointwise  estimate  exceeds  the 
tolerance  "tol",  finer  grids  are  added  as  leaf  nodes  to  the  tree.  This 
procedure  is  similar  to  one  used  by  Berger  (ref  4);  however,  there  are  more 
economical  error  estimation  strategies  (see,  for  example,  Bieterman  and 
Babuska  (refs  6,7))  which  we  are  currently  investigating. 

In  order  to  solve  the  finite  element  system,  Eq.  (3),  we  need  to  supply 
initial  and  boundary  conditions.  On  any  grid  with  p  ■  0,  a  =  a,  or  3  =  b, 
these  can  be  obtained  from  the  Initial  condition,  Eq.  (2),  or  prescribed 
boundary  conditions.  However,  artificial  initial  and  boundary  conditions  must 
be  created  at  all  other  coarse-fine  mesh  interfaces.  This  is  a  difficult  and 
crucial  problem  that  is  discussed  for  explicit  finite  difference  methods  by 

^Berger,  M.  J.,  "Adaptive  Mesh  Refinement  for  Hyperbolic  Partial  Differential 
Equations,"  Report  No.  STAN-CS-82-924,  Department  of  Computer  Science, 
Stanford  University,  1982. 

^Bieterman,  M.  and  Babuska,  I.,  "The  Finite  Element  Method  for  Parabolic 
Equations,  I.  A  Posteriori  Estimation,"  Numer.  Math.,  Vol.  40,  1982,  pp. 
339-371. 

^Bieterman,  M.  and  Babuska,  I. ,  "The  Finite  Element  Method  for  Parabolic 
Equations,  II.  A  Posteriori  Error  Estimation  and  Adaptive  Approach,"  Nume r . 
Math. .  Vol.  40,  1982,  pp.  373-406. 


Berger  (refs  4,8);  however,  it  is  largely  unanswered  for  finite  element 
applications.  Instabilities  or  incorrect  solutions  (see  Example  1  in  the 
following  section)  can  result  if  inappropriate  conditions  are  specified. 

For  initial  conditions,  two  strategies  immediately  come  to  mind:  (1) 
saving  all  fine  grid  data  for  propagation  in  time,  or  (2)  interpolating  the 
best  coarse  grid  data  to  finer  grids.  We  consider  a  blend  of  the  two 
strategies  which  consists  of  saving  the  fine  grid  data  down  to  a  given  level  X 
in  the  tree  and  subsequently  interpolating  for  finer  grids.  Each  grid  in  the 
first  X  levels  either  has  a  linked  list  of  the  initial  data  directly 
associated  with  it  or  uses  an  initial  data  list  of  an  ancestor  grid.  To  find 
the  value  of  the  solution  at  some  new  initial  point,  the  coordinate  of  that 
point  is  sequentially  compared  to  values  in  the  linked  list  until  an  interval 
containing  the  point  is  found  so  that  interpolation  can  be  used.  This  is 
costly  and  we  are  investigating  more  efficient  procedures  that  use  the  natural 
ordering  that  already  exist.  We  used  either  piecewise  linear  interpolation 
or  piecewise  parabolic  interpolation  with  shape  preserving  splines  developed 
by  McLaughlin  (ref  9).  For  each  grid  in  the  first  X  levels  of  the  tree,  a 
linked  list  is  created  to  store  the  initial  data.  We  are  studying  several 
alternative  ways  of  determining  a  proper  value  for  X. 


^Berger,  M.  J.,  "Adaptive  Mesh  Refinement  for  Hyperbolic  Partial  Differential 
Equations,"  Report  No.  STAN-CS-82-924,  Department  of  Computer  Science, 
Stanford  University,  1982. 

®Berger,  M.  J.,  "Stability  of  Interfaces  with  Mesh  Refinement,"  Report  No. 
83-42,  Institute  for  Computer  Applications  in  Science  and  Engineering,  NASA 
Langley  Research  Center,  Hampton,  1983. 

^McLaughlin,  H.  W.,  "Shape  Preserving  Planar  Interpolation:  An  Algorithm," 
IEEE  Computer  Graphics  and  Applies.,  Vol.  3,  1983,  pp.  58-67. 
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At  the  present  time,  we  prescribe  internal  Dirichlet  boundary  conditions 
by  linearly  interpolating  from  coarse  to  finer  grids.  A  buffer  zone  of  two 
elements  is  added  to  each  end  of  regions  of  high  error  that  do  not  intersect 
the  boundaries  x  =*  a  and  b.  If  two  buffer  zones  overlap  or  are  separated  from 
one  another  by  one  element,  the  two  grids  are  joined.  Similarly,  if  the 
buffer  is  only  one  element  away  from  either  a  or  b,  that  element  is  added  to 
the  grid. 


NUMERICAL  EXAMPLES 


An  experimental  code  based  on  the  algorithms  in  the  previous  section  has 
been  written  in  FORTRAN-77.  We  are  testing  it  on  several  examples,  some  of 
these  follow  and  others  are  presented  in  Reference  3.  All  results  were 
computed  in  double  precision  on  an  IBM  3081D  computer. 

Example  1 

In  order  to  illustrate  the  importance  of  adequately  resolving  initial 
conditions  at  each  time  step,  we  solve  the  linear  hyperbolic  initial  value 
problem 


u(x,0)  ■  u°(x)  ■ 


ut  +  ux  -  0  , 

(1/2) (cos(20  tt(x-0.45)-1)  ,  0.35  <  x  <  0.75 


^0  ,  otherwise 

We  3olve  this  problem  for  one  coarse  time  step  of  At  *  0.05,  10  elements  on  0 
<  x  <  1,  tol  *  0.01.  For  small  enough  times,  the  exact  solution  is  u°(x-t). 


^Flaherty,  J.  E.  and  Moore,  P.  K. ,  "An  Adaptive  Local  Refinement  Finite 
Element  Method  for  Parabolic  Partial  Differential  Equations,"  Proceedings  of 
the  International  Conference  on  Accuracy  and  Estimates  and  Adaptive 
Refinements  in  Finite  Element  Computations,  Technical  University  of  Lisbon, 
Lisbon,  Vol.  2,  1984,  pp.  139-152. 
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If  Initial  conditions  are  interpolated  from  the  coarse  to  the  fine  grid,  the 
oscillations  are  missed  and  an  incorrect  solution  is  computed,  possibly 
without  a  user  realizing  that  there  is  anything  wrong.  However,  saving 
initial  values  for  the  first  eight  level?  of  the  tree  of  grids  calculates  the 
correct  solution  to  the  prescribed  accuracy.  The  incorrect  and  correct 
solutions  are  shown  at  t  *  0.05  in  Figure  2. 

Example  2 

We  solve  the  following  problem  for  Burgers'  equation: 

ut  +  uux  =*  duxx  ,0<x<l  ,  0  <  t  <  1 
u(x,0)  =  sin  me  ,  0  <  x  <  1 
u(0,t)  -  u( 1 , t )  -  0  ,  t  >  0 

We  choose  d  ■  0.00003,  a  coarse  grid  of  10  elements  and  At  ■  0.1,  and 
piecewise  parabolic  approximations  for  the  initial  conditions  with  X  -  6.  It 
is  well  known  that  the  solution  of  this  problem  is  a  "pulse"  that  steepens  as 
it  travels  to  the  right  until  it  forms  a  shock  layer  at  x  -  1.  After  a  time 
of  0(l/d)  the  pulse  dissipates  and  the  solution  decays  to  zero.  We  solve  this 
problem  for  tol  ■  0.01  and  0.001  and  show  the  solutions  at  t  -  0.4  in  Figure 
3.  The  solution  with  the  cruder  tolerance  is  exhibiting  some  oscillations 
that  are  within  our  bounds.  These,  however,  are  not  visible  when  the  finer 
tolerance  is  used  to  solve  the  problem. 

Example  3 

We  solve  the  model  combustion  problem 

ut  +  ux  -  2eu  -  uxx  ,  0  <  x  <  1  ,  0  <  t  <  1 


u(x,0)  -  0  ,  u(0,t)  -  0  ,  UjjCl.t)  -  0 


The  exponential  nonlinearity  is  typical  in  combustion  problems  having 
Arrhenius  chemical  kinetics.  However,  in  this  case  the  solution  develops  a 
"hot  spot"  at  x  ■  1  and  becomes  infinite  when  t  is  approximately  0.85.  We 
choose  a  coarse  grid  of  20  elements  and  At  *  0.05,  tol  ■  0.001,  and  X  ■  6.  In 
Figure  4  we  show  the  computed  solution  U(x,t)  as  a  function  of  x  for  t  ■  0.05, 
0.06,  and  0.08,  and  in  Figure  5  we  show  the  mesh  that  was  used  to  solve  the 
problem.  We  see  that  the  mesh  is  initially  concentrated  in  the  region  near  x 
=*  0  where  the  curvature  of  the  solution  is  largest.  As  time  progresses  and 
the  curvature  diminishes,  excessive  refinement  is  not  necessary.  Finally,  as 
the  solution  begins  to  "blow-up”  our  algorithm  generates  a  fine  mesh  only  in 
the  region  near  x  =  1. 

DISCUSSIONS  AND  CONCLUSIONS 

We  have  briefly  described  an  adaptive  local  refinement  algorithm  for 
solving  time  dependent  partial  differential  equations.  Even  though  this  is 
very  much  a  working  algorithm  and  not  a  production  code,  we  are  very 
encouraged  by  the  preliminary  results.  We  are  investigating  several  possible 
ways  of  improving  the  efficiency  and  robustness  of  our  algorithm.  These 
include  adding  higher  order  polynomial  finite  element  approximations, 
adaptively  changing  the  number  of  elements  that  are  carried  forward  in  the 
coarse  grid  at  each  coarse  time  step,  selecting  the  appropriate  buffer 
length,  adaptively  determining  the  optimal  number  of  levels  of  initial 
conditions  at  coarse-fine  interfaces,  and  applying  the  best  boundary 
conditions  to  apply  at  internal  boundaries.  We  are  encouraged  by  the 
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performance  of  McLaughlin's  (ref  9)  shape  preserving  parabolic  splines; 
however,  the  entire  area  of  interpolating  from  coarse  to  fine  grids  needs 
further  study.  We  are  also  developing  non-Dirichlet  "natural"  boundary 
conditions  to  use  at  coarse-fine  mesh  interfaces. 

Finally,  we  are  very  interested  in  combining  the  moving  mesh  strategy  of, 
for  example,  References  5  and  10  with  the  present  local  refinement  strategy 
and  extending  our  methods  to  two  and  three  dimensions. 


^Davis,  S.  F.  and  Flaherty,  J.  E.,  "An  Adaptive  Finite  Element  for  Initial- 
Boundary  Value  Problems  for  Partial  Differential  Equations,"  SIAM  J.  Scl. 
Stat.  Comput. ,  Vol.  3,  1982,  pp.  6-27. 

^McLaughlin,  H.  W.,  "Shape  Preserving  Planar  Interpolation:  An  Algorithm," 
IEEE  Computer  Graphics  and  Applies.,  Vol.  3,  1983,  pp.  58-67. 
lyFlaherty,  J.  E.,  Coyle,  J.  M. ,  Ludwig,  R.,  and  Davis,  S.  F.,  "Adaptive 
Finite  Element  Methods  for  Parabolic  Partial  Differential  Equations,"  in 
Adaptive  Computational  Methods  for  Partial  Differential  Equations,  Babuska, 
I.,  Chandra,  J.,  and  Flaherty,  J.  E.  (eds.),  SIAM,  Philadelphia,  1983. 
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Figure  2.  Solution  of  Example  1  at  time  t  =  0.05  using  interpolation  from 
the  coarse  grid  to  the  fine  grid  (top)  and  saving  the  initial  values  for  the 
first  8  levels  of  the  tree  (bottom).  The  upper  solution  overlooks  the 
oscillations  and  is  incorrect. 
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